//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>

#include "ephi.hpp"

static FILE *fp;
typedef std::vector<vect3d> vlist_t;
static vlist_t planes;
static vlist_t points;
static std::vector<transf3d> tr;

static bool split (const vlist_t& vl, const vect3d& n, vlist_t& l1, vlist_t& l2)
{
  l1.clear();
  l2.clear();
  for ( size_t i = 0; i < vl.size(); i++ )
    {
      if ( n * vl[i] >= 0 )
          l1.push_back (vl[i]);
      else
          l2.push_back (vl[i]);
    }
  return !l1.empty() && !l2.empty();
}

static char *space (size_t n)
{
  static char buf[128];
  assert (n < sizeof(buf) - 1);
  memset (buf, ' ', n);
  buf[n] = 0;
  return buf;
}

static const char *fmt (prec_t f)
{
  static char buf[128];
  if ( fabs(fabs(f - prec_t_trunc(f)) - PREC_SQRT2 + 1) < 0.0001 )
    {
      prec_t m = prec_t_trunc(f);
      if ( m < 0 )
          m += 1;
      else
          m -= 1;
      if ( fabs (m) < 0.0001 )
          sprintf (buf, "%sPERC_SQRT2", f < 0 ? "-" : "");
      else
          sprintf (buf, "%f %s PERC_SQRT2", m, f < 0 ? "-" : "+");
    }
  else if ( fabs(fabs(f - prec_t_trunc(f)) - PREC_SQRT2 / 2) < 0.0001 )
    {
      if ( fabs(prec_t_trunc(f)) < 0.0001 )
          sprintf (buf, "%sPREC_SQRT2 / 2.0", f < 0 ? "-" : "");
      else
          sprintf (buf, "%f %s PREC_SQRT2 / 2.0", prec_t_trunc(f), f < 0 ? "-" : "+");
    }
  else if ( fabs(fabs(f) - 1) < 0.0001 )
    {
      sprintf (buf, "%s1", f < 0 ? "-" : "");
    }
  else if ( fabs(f) < 0.0001 )
    {
      sprintf (buf, "0");
    }
  else
    {
      sprintf (buf, "%f", f);
    }
  return buf;
}

static const char *format (const vect3d& v)
{
  static char buf[128];
  std::string r;
  prec_t a[3] = {v.x, v.y, v.z};
  for ( size_t i = 0; i < 3; i++ )
    {
      r += fmt(a[i]);
      if ( i != 2 ) r += ", ";
    }
  memcpy (buf, r.c_str(), r.size() + 1);
  return buf;
}

static const char *one (prec_t v)
{
  static char buf[128];
  assert (fabs(fabs(v)-1) < 0.0001);
  if ( v < 0 )
      strcpy (buf, "-");
  else
      buf[0] = 0;
  return buf;
}

static void compare_plane (const vect3d& plane)
{
  bool has = false;
  if ( fabs(plane.x) == fabs(plane.y) && plane.z == 0 ||
      fabs(plane.x) == fabs(plane.z) && plane.y == 0 ||
      fabs(plane.y) == fabs(plane.z) && plane.x == 0 )
    {
      if ( fabs(plane.x) > 0.1 )
        {
          fprintf (fp, "%sp.x", plane.x < 0 ? "-" : "");
          has = true;
        }

      if ( fabs(plane.y) > 0.1 )
        {
          if ( plane.y < 0 )
            {
              if ( has )
                  fprintf (fp, " - p.y");
              else
                  fprintf (fp, "-p.y");
            }
          else if ( has )
            {
              fprintf (fp, " + p.y");
            }
          else
            {
              fprintf (fp, "p.y");
            }
          has = true;
        }

      if ( fabs(plane.z) > 0.1 )
        {
          if ( plane.z < 0 )
            {
              if ( has )
                  fprintf (fp, " - p.z");
              else
                  fprintf (fp, "-p.z");
            }
          else if ( has )
            {
              fprintf (fp, " + p.z");
            }
          else
            {
              fprintf (fp, "p.z");
            }
        }
    }
  else
    {
      fprintf (fp, "p * vect3d(%s)\n", format(plane));
    }
  fprintf (fp, " >= 0");
}

static void recurse48 (const vlist_t& vl, size_t level)
{
  assert (vl.size() != 0);
  vlist_t v1, v2;
  if ( level == 0 )
    {
      split (vl, vect3d(1, 0, 0), v1, v2);
      fprintf (fp, "  if ( xge0 )\n");
      fprintf (fp, "    {\n");
      recurse48 (v1, 1);
      fprintf (fp, "    }\n");
      fprintf (fp, "  else\n");
      fprintf (fp, "    {\n");
      recurse48 (v2, 1);
      fprintf (fp, "    }\n");
    }
  else if ( level == 1 )
    {
      split (vl, vect3d(0, 1, 0), v1, v2);
      fprintf (fp, "      if ( yge0 )\n");
      fprintf (fp, "        {\n");
      recurse48 (v1, 2);
      fprintf (fp, "        }\n");
      fprintf (fp, "      else\n");
      fprintf (fp, "        {\n");
      recurse48 (v2, 2);
      fprintf (fp, "        }\n");
    }
  else if ( level == 2 )
    {
      split (vl, vect3d(0, 0, 1), v1, v2);
      fprintf (fp, "          if ( zge0 )\n");
      fprintf (fp, "            {\n");
      recurse48 (v1, 3);
      fprintf (fp, "            }\n");
      fprintf (fp, "          else\n");
      fprintf (fp, "            {\n");
      recurse48 (v2, 3);
      fprintf (fp, "            }\n");
    }
  else
    {
      // only one point left, but this entire sector can be tansformed
      // using this point's transformation matrix we worked out before
      if ( vl.size() == 1 )
        {
          for ( size_t i = 0; i < points.size(); i++ )
            {
              if ( vl[0] == points[i] )
                {
                  fprintf (fp, "%sp = vect3d(", space(2 + 4*level));
                  if ( fabs(tr[i].a) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].a));
                  if ( fabs(tr[i].b) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].b));
                  if ( fabs(tr[i].c) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].c));
                  fprintf (fp, ", ");
                  if ( fabs(tr[i].d) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].d));
                  if ( fabs(tr[i].e) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].e));
                  if ( fabs(tr[i].f) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].f));
                  fprintf (fp, ", ");
                  if ( fabs(tr[i].g) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].g));
                  if ( fabs(tr[i].h) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].h));
                  if ( fabs(tr[i].i) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].i));
                  fprintf (fp, ");\n");
                  fprintf (fp, "%sreturn transforms48[%u];\n", space(2 + 4*level), (unsigned int)i);
                  return;
                }
            }
          assert(0);
        }
      size_t i;
      for ( i = 0; i < planes.size(); i++ )
          if ( split (vl, planes[i], v1, v2) )
              break;
      assert (i != planes.size());
      if ( fabs(planes[i].x) == fabs(planes[i].y) && planes[i].z == 0 ||
          fabs(planes[i].x) == fabs(planes[i].z) && planes[i].y == 0 ||
          fabs(planes[i].y) == fabs(planes[i].z) && planes[i].x == 0 )
        {
          fprintf (fp, "%sif ( ", space(2 + 4*level));
          bool has = false;
          if ( fabs(planes[i].x) > 0.1 )
            {
              fprintf (fp, "%sp.x", planes[i].x < 0 ? "-" : "");
              has = true;
            }

          if ( fabs(planes[i].y) > 0.1 )
            {
              if ( planes[i].y < 0 )
                {
                  if ( has )
                      fprintf (fp, " - p.y");
                  else
                      fprintf (fp, "-p.y");
                }
              else if ( has )
                {
                  fprintf (fp, " + p.y");
                }
              else
                {
                  fprintf (fp, "p.y");
                }
              has = true;
            }

          if ( fabs(planes[i].z) > 0.1 )
            {
              if ( planes[i].z < 0 )
                {
                  if ( has )
                      fprintf (fp, " - p.z");
                  else
                      fprintf (fp, "-p.z");
                }
              else if ( has )
                {
                  fprintf (fp, " + p.z");
                }
              else
                {
                  fprintf (fp, "p.z");
                }
            }
          fprintf (fp, " >= 0 )\n");
        }
      else
        {
          fprintf (fp, "%sif ( p * vect3d(%s) >= 0)\n", space(2 + 4*level), format(planes[i]));
        }
      fprintf (fp, "%s  {\n", space(2 + 4*level));
      recurse48 (v1, level+1);
      fprintf (fp, "%s  }\n", space(2 + 4*level));
      fprintf (fp, "%selse\n", space(2 + 4*level));
      fprintf (fp, "%s  {\n", space(2 + 4*level));
      recurse48 (v2, level+1);
      fprintf (fp, "%s  }\n", space(2 + 4*level));
    }
}

static void recurse24 (const vlist_t& vl, size_t level)
{
  assert (vl.size() != 0);
  vlist_t v1, v2;
  if ( level == 0 )
    {
      split (vl, vect3d(1, 0, 0), v1, v2);
      fprintf (fp, "  if ( xge0 )\n");
      fprintf (fp, "    {\n");
      recurse24 (v1, 1);
      fprintf (fp, "    }\n");
      fprintf (fp, "  else\n");
      fprintf (fp, "    {\n");
      recurse24 (v2, 1);
      fprintf (fp, "    }\n");
    }
  else if ( level == 1 )
    {
      split (vl, vect3d(0, 1, 0), v1, v2);
      fprintf (fp, "      if ( yge0 )\n");
      fprintf (fp, "        {\n");
      recurse24 (v1, 2);
      fprintf (fp, "        }\n");
      fprintf (fp, "      else\n");
      fprintf (fp, "        {\n");
      recurse24 (v2, 2);
      fprintf (fp, "        }\n");
    }
  else if ( level == 2 )
    {
      split (vl, vect3d(0, 0, 1), v1, v2);
      fprintf (fp, "          if ( zge0 )\n");
      fprintf (fp, "            {\n");
      recurse24 (v1, 3);
      fprintf (fp, "            }\n");
      fprintf (fp, "          else\n");
      fprintf (fp, "            {\n");
      recurse24 (v2, 3);
      fprintf (fp, "            }\n");
    }
  else
    {
      if ( vl.size() == 1 )
        {
          for ( size_t i = 0; i < points.size(); i++ )
            {
              if ( vl[0] == points[i] )
                {
                  fprintf (fp, "%sp = vect3d(", space(2 + 4*level));
                  if ( fabs(tr[i].a) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].a));
                  if ( fabs(tr[i].b) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].b));
                  if ( fabs(tr[i].c) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].c));
                  fprintf (fp, ", ");
                  if ( fabs(tr[i].d) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].d));
                  if ( fabs(tr[i].e) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].e));
                  if ( fabs(tr[i].f) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].f));
                  fprintf (fp, ", ");
                  if ( fabs(tr[i].g) > 0.001 ) fprintf (fp, "%sp.x", one(tr[i].g));
                  if ( fabs(tr[i].h) > 0.001 ) fprintf (fp, "%sp.y", one(tr[i].h));
                  if ( fabs(tr[i].i) > 0.001 ) fprintf (fp, "%sp.z", one(tr[i].i));
                  fprintf (fp, ");\n");
                  fprintf (fp, "%sreturn transforms24[%u];\n", space(2 + 4*level), (unsigned int)i);
                  return;
                }
            }
          assert(0);
        }

      if ( vl.size() == 3 )
        {
          size_t i;
          for ( i = 0; i < planes.size(); i++ )
            {
              bool t1 = planes[i] * vl[0] >= 0.001;
              bool t2 = planes[i] * vl[1] >= 0.001;
              bool t3 = planes[i] * vl[2] >= 0.001;
              if ( t1 && !t2 && !t3 || !t1 && t2 && !t3 || !t1 && !t2 && t3 )
                  break;
            }
          assert (i != planes.size());
          vect3d plane = planes[i];
          for ( size_t i = 0; i < vl.size(); i++ )
            {
              if ( plane * vl[i] >= 0 || fabs(plane * vl[i]) <= 0.001 )
                  v1.push_back (vl[i]);
              if ( plane * vl[i] <= 0 || fabs(plane * vl[i]) <= 0.001 )
                  v2.push_back (vl[i]);
            }
          assert (v1.size() == 2 && v2.size() == 2);
          fprintf (fp, "%sif ( ", space(2+4*level));
          compare_plane (plane);
          fprintf (fp, " )\n");
          fprintf (fp, "%s  {\n", space(2+4*level));
          recurse24 (v1, level + 1);
          fprintf (fp, "%s  }\n", space(2+4*level));
          fprintf (fp, "%selse\n", space(2+4*level));
          fprintf (fp, "%s  {\n", space(2+4*level));
          recurse24 (v2, level + 1);
          fprintf (fp, "%s  }\n", space(2+4*level));
        }
      else
        {
          assert (vl.size() == 2 );
          size_t i;
          for ( i = 0; i < planes.size(); i++ )
            {
              if ( planes[i] * vl[0] >= 0.001 && planes[i] * vl[1] <= -0.001 ||
                  planes[i] * vl[0] <= -0.001 && planes[i] * vl[1] >= 0.001 )
                  break;
            }
          assert (i != planes.size());
          vect3d plane = planes[i];
          split (vl, plane, v1, v2);
          assert (v1.size() == 1 && v2.size() == 1);

          fprintf (fp, "%sif ( ", space(2+4*level));
          compare_plane (plane);
          fprintf (fp, " )\n");
          fprintf (fp, "%s  {\n", space(2+4*level));
          recurse24 (v1, level + 1);
          fprintf (fp, "%s  }\n", space(2+4*level));
          fprintf (fp, "%selse\n", space(2+4*level));
          fprintf (fp, "%s  {\n", space(2+4*level));
          recurse24 (v2, level + 1);
          fprintf (fp, "%s  }\n", space(2+4*level));
        }


#if 0
      // 3 points left. which is it?!?
      assert (vl.size() == 3);
      size_t i;
      for ( i = 0; i < planes.size(); i++ )
        {
          bool t1 = planes[i] * vl[0] >= 0.001;
          bool t2 = planes[i] * vl[1] >= 0.001;
          bool t3 = planes[i] * vl[2] >= 0.001;
          if ( t1 && !t2 && !t3 || !t1 && t2 && !t3 || !t1 && !t2 && t3 )
              break;
        }
      assert (i != planes.size());
      vect3d split = planes[i];

      const vect3d *cp = 0;
      const vect3d *lp = 0;
      const vect3d *p3 = 0;

      for ( size_t i = 0; i < 3; i++ )
        {
          prec_t r = vl[i] * split;
          printf ("%f\n", r);
          if ( r >= 0.001 )
              cp = &vl[i];
          else if ( fabs(r) <= 0.001 )
              lp = &vl[i];
          else
              p3 = &vl[i];
        }
      assert(cp != 0 && lp != 0 && p3 != 0);
      cp->print ("cp");

      for ( i = 0; i < planes.size(); i++ )
        {
          if ( *cp * planes[i] >= 0.001 && planes[i] != split )
              break;
        }
      assert (i != planes.size());
      vect3d split2 = planes[i];

      vlist_t next;

      fprintf (fp, "%sbool split1 = ", space(2+4*level));
      compare_plane (split);
      fprintf (fp, ";\n");
      fprintf (fp, "%sbool split2 = ", space(2+4*level));
      compare_plane (split2);
      fprintf (fp, ";\n");
      fprintf (fp, "%sif ( split1 && split2 )\n", space(2+4*level));
      fprintf (fp, "%s  {\n", space(2+4*level));
      next.clear();
      next.push_back (*cp);
      recurse24 (next, level + 1);
      // we have cp index
      fprintf (fp, "%s  }\n", space(2+4*level));
      fprintf (fp, "%selse if (split1 && !split2)\n", space(2+4*level));
      fprintf (fp, "%s  {\n", space(2+4*level));
      next.clear();
      next.push_back (*lp);
      recurse24 (next, level + 1);
      fprintf (fp, "%s  }\n", space(2+4*level));
      fprintf (fp, "%selse if (split2)\n", space(2+4*level));
      fprintf (fp, "%s  {\n", space(2+4*level));
      next.clear();
      next.push_back (*p3);
      recurse24 (next, level + 1);
      fprintf (fp, "%s  }\n", space(2+4*level));

      printf ("Boo!\n");
#endif
    }
}

static void make24 ()
{
  planes.clear();
  planes.push_back (vect3d(1, 1, 0).normal());
  planes.push_back (vect3d(1, -1, 0).normal());
  planes.push_back (vect3d(1, 0, 1).normal());
  planes.push_back (vect3d(1, 0, -1).normal());
  planes.push_back (vect3d(0, 1, 1).normal());
  planes.push_back (vect3d(0, 1, -1).normal());
  planes.push_back (-vect3d(1, 1, 0).normal());
  planes.push_back (-vect3d(1, -1, 0).normal());
  planes.push_back (-vect3d(1, 0, 1).normal());
  planes.push_back (-vect3d(1, 0, -1).normal());
  planes.push_back (-vect3d(0, 1, 1).normal());
  planes.push_back (-vect3d(0, 1, -1).normal());

  tr.clear();
  tr.push_back (transf3d ( 1,  0,  0,  0,  1,  0,  0,  0,  1));
  tr.push_back (transf3d ( 1,  0,  0,  0, -1,  0,  0,  0, -1));
  tr.push_back (transf3d ( 1,  0,  0,  0,  0, -1,  0,  1,  0));
  tr.push_back (transf3d ( 1,  0,  0,  0,  0,  1,  0, -1,  0));
  tr.push_back (transf3d (-1,  0,  0,  0,  1,  0,  0,  0, -1));
  tr.push_back (transf3d (-1,  0,  0,  0, -1,  0,  0,  0,  1));
  tr.push_back (transf3d (-1,  0,  0,  0,  0,  1,  0,  1,  0));
  tr.push_back (transf3d (-1,  0,  0,  0,  0, -1,  0, -1,  0));

  tr.push_back (transf3d ( 0,  1,  0,  0,  0,  1,  1,  0,  0));
  tr.push_back (transf3d ( 0, -1,  0,  0,  0, -1,  1,  0,  0));
  tr.push_back (transf3d ( 0,  0, -1,  0,  1,  0,  1,  0,  0));
  tr.push_back (transf3d ( 0,  0,  1,  0, -1,  0,  1,  0,  0));
  tr.push_back (transf3d ( 0,  1,  0,  0,  0, -1, -1,  0,  0));
  tr.push_back (transf3d ( 0, -1,  0,  0,  0,  1, -1,  0,  0));
  tr.push_back (transf3d ( 0,  0,  1,  0,  1,  0, -1,  0,  0));
  tr.push_back (transf3d ( 0,  0, -1,  0, -1,  0, -1,  0,  0));

  tr.push_back (transf3d ( 0,  0,  1,  1,  0,  0,  0,  1,  0));
  tr.push_back (transf3d ( 0,  0, -1,  1,  0,  0,  0, -1,  0));
  tr.push_back (transf3d ( 0,  1,  0,  1,  0,  0,  0,  0, -1));
  tr.push_back (transf3d ( 0, -1,  0,  1,  0,  0,  0,  0,  1));
  tr.push_back (transf3d ( 0,  0, -1, -1,  0,  0,  0,  1,  0));
  tr.push_back (transf3d ( 0,  0,  1, -1,  0,  0,  0, -1,  0));
  tr.push_back (transf3d ( 0,  1,  0, -1,  0,  0,  0,  0,  1));
  tr.push_back (transf3d ( 0, -1,  0, -1,  0,  0,  0,  0, -1));

  fprintf (fp, "static transf3d_pair transforms24[24] = {\n");
  for ( size_t i = 0; i < tr.size(); i++ )
    {
      transf3d& t = tr[i];
      transf3d inv = t.inverse();
      fprintf (fp, "  {transf3d(");
      fprintf (fp, "%s, ", fmt(t.a));
      fprintf (fp, "%s, ", fmt(t.b));
      fprintf (fp, "%s, ", fmt(t.c));
      fprintf (fp, "%s, ", fmt(t.d));
      fprintf (fp, "%s, ", fmt(t.e));
      fprintf (fp, "%s, ", fmt(t.f));
      fprintf (fp, "%s, ", fmt(t.g));
      fprintf (fp, "%s, ", fmt(t.h));
      fprintf (fp, "%s", fmt(t.i));
      fprintf (fp, "), transf3d(");

      fprintf (fp, "%s, ", fmt(inv.a));
      fprintf (fp, "%s, ", fmt(inv.b));
      fprintf (fp, "%s, ", fmt(inv.c));
      fprintf (fp, "%s, ", fmt(inv.d));
      fprintf (fp, "%s, ", fmt(inv.e));
      fprintf (fp, "%s, ", fmt(inv.f));
      fprintf (fp, "%s, ", fmt(inv.g));
      fprintf (fp, "%s, ", fmt(inv.h));
      fprintf (fp, "%s", fmt(inv.i));

      fprintf (fp, ")}%s\n", i != tr.size() - 1 ? "," : "");
    }
  fprintf (fp, "};\n");
  fprintf (fp, "\n");

  points.clear();
  for ( size_t i = 0; i < tr.size(); i++ )
      points.push_back (tr[i].inverse() * vect3d (1, 0.1, 0.1));

  fprintf (fp, "const transf3d_pair& cube_polywell_transform24 (vect3d& p)\n");
  fprintf (fp, "{\n");
  fprintf (fp, "  bool xge0 = p.x >= 0;\n");
  fprintf (fp, "  bool yge0 = p.y >= 0;\n");
  fprintf (fp, "  bool zge0 = p.z >= 0;\n");

  recurse24 (points, 0);

  fprintf (fp, "}\n");
  fprintf (fp, "\n");
  fprintf (fp, "const transf3d_pair* get_cube_polywell_transforms24 ()\n");
  fprintf (fp, "{\n");
  fprintf (fp, "  return transforms24;\n");
  fprintf (fp, "}\n");
  fprintf (fp, "\n");

#if 0
  vlist_t ml;
  ml.push_back (vect3d(0, 1, 0));
  ml.push_back (vect3d(0, 0, 1));
  ml.push_back (vect3d (1, -1, 0).normal());
  ml.push_back (vect3d (1, 0, -1).normal());

  for ( size_t i = 0; i < points.size(); i++ )
    {
      size_t j;
      for ( j = 0; j < ml.size(); j++ )
          if ( ml[j] * points[i] < 0.001 )
              break;
      if ( j != ml.size() )
          continue;
      printf ("match: %u\n", (unsigned int)i);
    }

  points[0].print("points[0]");
  for ( size_t i = 0; i < planes.size(); i++ )
    {
      prec_t f = planes[i] * points[0];
      char b[128];
      if ( f > 0.001 )
          printf ("%f: %s\n", f, planes[i].sprint(b));
    }
#endif
}

int main (int argc, char *argv[])
{
  // now I wish I had paid attention at the linear algebra class,
  // but I did not, so I have to resort to this

  // there are 3 + 6 splitting planes
  planes.push_back (vect3d(1, 1, 0).normal());
  planes.push_back (vect3d(1, -1, 0).normal());
  planes.push_back (vect3d(1, 0, 1).normal());
  planes.push_back (vect3d(1, 0, -1).normal());
  planes.push_back (vect3d(0, 1, 1).normal());
  planes.push_back (vect3d(0, 1, -1).normal());

  // there are 48 symmetric sectors, generate 48 points within those sectors

  points.clear();
  // 8 points in each face
  for ( size_t i = 0; i < 8; i++ )
    {
      prec_t phi = 1 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (3, a, b));
    }

  for ( size_t i = 0; i < 8; i++ )
    {
      prec_t phi = 3 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (-3, a, b));
    }

  for ( size_t i = 0; i < 8; i++ )
    {
      prec_t phi = 1 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (a, b, 3));
    }

  for ( size_t i = 0; i < 8; i++ )
    {
      prec_t phi = 3 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (a, b, -3));
    }

  for ( size_t i = 0; i < 8; i++ ) // correct
    {
      prec_t phi = 3 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (a, 3, b));
    }

  for ( size_t i = 0; i < 8; i++ )
    {
      prec_t phi = 1 * (2 * PREC_PI / 16) + i * 2 * PREC_PI / 8.0;
      prec_t a = cos(phi);
      prec_t b = sin(phi);
      points.push_back (vect3d (a, -3, b));
    }

  fp = fopen ("polywell_sectors.cpp", "w");
  if ( !fp )
      exit(-1);

  fprintf (fp, "// This file was automatically generated by sectors.cpp\n");
  fprintf (fp, "// This file is part of ephi (http://www.mare.ee/indrek/ephi/)\n");
  fprintf (fp, "\n");
  fprintf (fp, "#include \"math3d.hpp\"\n");
  fprintf (fp, "\n");
  fprintf (fp, "// 3x3 transformation matrixes, forward and backwards\n");
  fprintf (fp, "static transf3d_pair transforms48[48] = {\n");

  vlist_t limits;
  limits.push_back (vect3d(1, 0, 0));
  limits.push_back (vect3d(0, 1, 0));
  limits.push_back (vect3d(0, 0, 1));
  limits.push_back (vect3d(1, 0, -1).normal());
  limits.push_back (vect3d(0, -1, 1).normal());

  for ( size_t i = 0; i < points.size(); i++ )
    {
      transf3d t;
      t.identity();
      vect3d p = points[i];
      while ( true )
        {
          size_t j;
          for ( j = 0; j < limits.size(); j++ )
            {
              if ( p * limits[j] < 0 )
                  break;
            }

          if ( j == limits.size() )
              break;

          // reflect point on the other side and merge final transformation matrix
          transf3d r;
          r.reflect (limits[j]);
          //vect3d oldp = p;
          p = r * p;
          //char b1[128], b2[128], b3[128];
          //printf ("%u: reflect on %s, %s->%s\n", (unsigned int)i, limits[j].sprint(b1), oldp.sprint(b2), p.sprint(b3));
          //t.print("t");
          //r.print("r");
          t = r * t;
          //t.print("result");
        }

      vect3d tp = t * points[i];

      transf3d inv = t.inverse();
      vect3d ip = inv * tp;

#if 0
      points[i].print ("point");
      p.print ("mp");
      tp.print ("tp");
      ip.print ("ip");

      t.print ("  t");
      printf ("t.det=%f\b", t.det());
      inv.print ("inv");
      (inv * t).print("inv * t");
      (t * inv).print("t * inv");
#endif

      assert ( (tp - p).length() < 0.0001); // checking the transformation works
      assert ( (points[i] - ip).length() < 0.0001); // checking the inversion works as well

      tr.push_back (t);

      fprintf (fp, "  {transf3d(");
      fprintf (fp, "%s, ", fmt(t.a));
      fprintf (fp, "%s, ", fmt(t.b));
      fprintf (fp, "%s, ", fmt(t.c));
      fprintf (fp, "%s, ", fmt(t.d));
      fprintf (fp, "%s, ", fmt(t.e));
      fprintf (fp, "%s, ", fmt(t.f));
      fprintf (fp, "%s, ", fmt(t.g));
      fprintf (fp, "%s, ", fmt(t.h));
      fprintf (fp, "%s", fmt(t.i));
      fprintf (fp, "), transf3d(");

      fprintf (fp, "%s, ", fmt(inv.a));
      fprintf (fp, "%s, ", fmt(inv.b));
      fprintf (fp, "%s, ", fmt(inv.c));
      fprintf (fp, "%s, ", fmt(inv.d));
      fprintf (fp, "%s, ", fmt(inv.e));
      fprintf (fp, "%s, ", fmt(inv.f));
      fprintf (fp, "%s, ", fmt(inv.g));
      fprintf (fp, "%s, ", fmt(inv.h));
      fprintf (fp, "%s", fmt(inv.i));

      fprintf (fp, ")}%s\n", i != points.size() - 1 ? "," : "");
    }

  fprintf (fp, "};\n");
  fprintf (fp, "\n");

  fprintf (fp, "// translates a generic point into the specific sector in the 1/48 space\n");
  fprintf (fp, "// the specific sector is limited by following planes:\n");
  fprintf (fp, "// (1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 0, -1), (0, -1, 1)\n");
  fprintf (fp, "// returns the transformation matrix used and also an inverse of the matrix\n");
  fprintf (fp, "// that can be used to translate back\n");

  fprintf (fp, "const transf3d_pair& cube_polywell_transform48 (vect3d& p)\n");
  fprintf (fp, "{\n");
  fprintf (fp, "  bool xge0 = p.x >= 0;\n");
  fprintf (fp, "  bool yge0 = p.y >= 0;\n");
  fprintf (fp, "  bool zge0 = p.z >= 0;\n");

  recurse48 (points, 0);

  fprintf (fp, "}\n");
  fprintf (fp, "\n");
  fprintf (fp, "const transf3d_pair* get_cube_polywell_transforms48 ()\n");
  fprintf (fp, "{\n");
  fprintf (fp, "  return transforms48;\n");
  fprintf (fp, "}\n");
  fprintf (fp, "\n");

  make24();
  fclose (fp);

  return 0;
}

